Arching in tapped deposits of hard disks. 
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We simulate the tapping of a bed of hard disks in a rectangular box by using a pseudo dynamic 
algorithm. In these simulations, arches are unambiguously denned and we can analyze their proper- 
ties as a function of the tapping amplitud. We find that an order-disorder transition occurs within 
a narrow range of tapping amplitudes as it has been seen by others. Arches are always present in 
the system althought they exhibit regular shapes in the ordered regime. Interestingly, an increase 
in the number of arches does not always correspond to a reduction in the packing fraction. This is 
in contrast with what is found in three-dimensional systems. 



I. INTRODUCTION 

The study of structural properties of assemblies of 
macroscopic particles packed in a container has become 
an important area of granular matter. Properties such as 
packing fraction (</>), coordination number (z), and arch 
distributions depend on the way the packings are created. 
The so-called " Chicago experiment" pj has shown that a 
tapped granular bed can achieve a stationary state where 
the system characterizing parameters (in this case <j>) de- 
pend only on the tapping amplitude. This simplifies our 
investigations when we are interested only in the steady 
state of the system. 

Arching is one of the collective phenomena that many 
associate to the appearance of voids in a granular sam- 
ple that leads to the lowering of <j>. Moreover, arching is 
directly related to the reduction of particle-particle con- 
tacts in the assembly, which determines the value of z 
[EIE|< Arch formation is crucial in the jamming of gran- 
ular flows 0, IE IE an d ^ nas a l so b een proposed as a 
mechanism for size segregation [E @ ■ 

Two-dimensional (2D) granular systems merit atten- 
tion due to the multiple applications in industry and ur- 
ban life (pills, bottles, etc. on a conveyor belt [10||, traffic 
jams ^3 > an d crowd control 01) an d a lso due to the pos- 
sibility of testing theoretical models that are particularly 
simple to treat in 2D. 

Simulation studies of the shaking of 2D granular beds 
have been carried out by others in the past (see for ex- 
ample Ref. 0). However, these pseudo dynamic sim- 
ulations do not consider simultaneous deposition of the 
grains. This prevents arching during deposition, which 
is a main ingredient in granular systems. On the other 
hand, simulations through realistic granular dynamics 
have the drawback that there is not a clear criterion to 
decide when the system is fully deposited, since kinetic 
energy does not vanishes completely in these simulations. 
Besides, identifying which particles support a given par- 
ticle in these type of simulations may prove rather com- 
plex. 

In this work we study the formation of arches in a 
2D assembly of hard disks by using a pseudo dynamic 
simulation of the tapping and simultaneous deposition of 



inelastic rough particles. We follow the changes in 0, z, 
and arch properties as a function of the tapping ampli- 
tude. Also, an annealed tapping similar to the "Chicago 
experiment" is carried out on the system. 



II. ARCHING AND COORDINATION NUMBER 

Arches are multi-particle structures formed during si- 
multaneous (non-sequential) deposition of an assembly 
of granular particles. Particles in an arch support each 
other so that the whole set remains stable againstgravity 
(or the driving force that promotes deposition) 

In a two-dimensional bed of convex particles con- 
structed by sequential deposition, the mean coordination 
number (z) is 4: each newly deposited particle adds two 
contacts — which stabilize the particle — to the system. 
If particles are deposited non-sequentially, mutual con- 
tacts are created where two particles share a stabilizing 
contact. This diminishes in 1 — with respect to a se- 
quential deposition — the total number of contacts in the 
system and so the value of (z) drops. Arches in 2D are 
string-like: any arch has two end-particles that share a 
single mutual contact and a variable number of center- 
particles that share two mutual contacts each. If n(s) is 
the number of arches consisting of s particles in a system 
of N particles, there are 
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center-particles in the assembly. Here, n(s — 1) corre- 
sponds to the number of particles that do not belong to 
any arch. Therefore, the mean coordination number can 
be obtained as 
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We use the subscript "support" because this coordina- 
tion number does not take into account the existence of 
non-supporting contacts. In non-sequentially deposited 
beds there will be some contacts that do not serve to the 
stability of any of the two touching particles. In our sim- 
ulations, these non-supporting contacts represent least 
than 0.3% of the contacts. Using Eqs. i an< 4 after 

some algebra we obtain 
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The analysis above can be done also in 3D but it re- 
quires more complex expressions due to the fact that 
arches can present branches 0, 0, 0] • 

Bearing in mind that X^fLi sn(s) = N, we can model 
n(s) as a simple exponential decay in the limit N — > oo, 
i.e., 



n(s)/N = 2(cosh[a] - l)exp[-as], (5) 
with a a positive number. In this model we obtain 
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The exponential model seems to work well for string- 
like bridges in 3D (seeRef. Q). We will see below that, 
in 2D, this is only a rough model for n(s) according to our 
simulation results. However, from Eq. 0we can already 
observe that (z) support can only vary between 2 and 4 in 
2D. Notice that this result does not depend on the shape 
of the particles so long as they are convex and support 
each other through point contacts. The case (z) suppor t = 
2 corresponds to a very particular configuration where all 
particles belong to a single large arch. This, of course, is 
impossible to achieve with walls present. 

In a graph of p-z/N as a function of pi/N (see 
Fig. 1) any packing must belong to the triangle 
[(0,0), (0, 1), (1,0)] and the lines parallel to the segment 
[(0, 0.5), (1, 0)] correspond to states of equal (z) S up P ort- 
The curve shown in Fig. 1 corresponds to the loci of 
the states given by the exponential model of n(s) with a 
ranging from to oo. In practice, (z) S upport only ranges 
from 3 to 4 (see the simulation data point in Fig. 1). 
The lower limit ((z) support = 3) has been suggested to 
correspond to the marginally rigid state in 2D 



III. SIMULATION DETAILS 

Our simulations are based on an algorithm for inelas- 
tic hard disks designed by Manna and Khakhar 0, . 
This is a pseudo dynamics that consists in small falls and 
rolls of the grains until they come to rest by contacting 
other particles or the system boundaries. We use a con- 
tainer formed by a flat base and two flat vertical walls. 



No periodic boundary conditions are applied. Once all 
the grains come to rest, the system is expanded and ran- 
domly shaken to simulate a vertical tap. Then, a new 
deposition cycle begins. After several taps, the system 
achieves a steady state where all characterizing parame- 
ters fluctuate around an equilibrium value. 

The deposition algorithm consists in picking up a disk 
in the system and performing a free fall of length S if the 
disk has no supporting contacts, or a roll of arc-length 5 
over its supporting disk if the disk has one single support- 
ing contact jl7t fl8| . Disks with two supporting contacts 
are considered stable and left in their positions. If in 
the course of a fall of length 5 a disk collides with an- 
other disk (or the base), the falling disk is put just in 
contact and this contact is defined as its first supporting 
contact. Analogously, If in the course of a roll of length 
5 a disk collides with another disk (or a wall), the rolling 
disk is put just in contact. If the first supporting con- 
tact and the second contact are such that the disk is in a 
stable position, the second contact is defined as the sec- 
ond supporting contact; otherwise, the lowest of the two 
contacting particle is taken as the first supporting con- 
tact of the rolling disk and the second supporting contact 
is left undefined. If, during a roll, a particle reaches a 
lower position than the supporting particle over which it 
is rolling, its first supporting contact is left undefined. A 
moving disk can change the stability state of other disks 
supported by it, therefore, this information is updated 
after each move. The deposition is over once each par- 
ticle in the system has both supporting contacts defined 
or is in contact with the base (particles at the base are 
supported by a single contact). Then, the coordinates of 
the centers of the disks and the corresponding labels of 
the two supporting particles, wall, or base, are saved for 
analysis. 

The tapping of the system is simulated by multiply- 
ing the vertical coordinate of each particle by A (with 
A > 1). Then, the particles are subjected to sev- 
eral (about 20) Monte Carlo loops where positions are 
changed by displacing particles a random length Ar uni- 
formly distributed in the range < Ar < A — 1. New 
configurations that correspond to overlaps are rejected. 
This disordering phase is crucial to avoid particles falling 
back again into the same positions. Moreover, the upper 
limit for Ar (i.e. A — 1) is deliberately chosen so that a 
larger tap promotes larger random changes in the particle 
positions. 

The simulations are carried out in a rectangular box 
of width 20 containing 2000 equal-sized disks of radius 
r = 0.1 + l/y/2 w 0.807. The disk size is chosen in such 
way that the neighbor linked-cell system used to speed up 
the simulation contains a single particle per cell. Initially, 
disks are placed at random in the simulation box (with 
no overlaps) and deposited using the pseudo dynamic 
algorithm. Then, 10 3 tapping cycles are performed for 
equilibration followed by 10 3 taps for production; saving 
only 1 every 10 fully deposited configurations. Tapping 
amplitudes range from 1.02 up to 2.0. Enlarging the 
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width of the simulation box up to 30 has shown no effect 
on the results that we present here. The same is true if 
one reduces the number of particles down to 500. 

An important point in this simulations is the effect 
that the parameter 8 has in the results since particles do 
not move simultaneously but one at a time. One might 
expect that in the limit 5 — ► we should recover a fairly 
" realistic" dynamics for fully inelastic rough disk dragged 
downwards at constant velocity. This should represent 
particles deposited in a viscous medium or carried by a 
conveyor belt. In Fig. 2, we show <ft and (z) SU pport for 
A = 1.1 as a function of S. As we can see, the results 
are independent of S for small values of the parameter. 
Since computation efficiency decreases with decreasing S 
— due to the number of free falls required for the particles 
to come together at the bottom of the container — , we 
choose the largest value that yield results indistinguish- 
able from the small-<5 limit within statistical errors (i.e., 
5 = 0.01). This value might be inappropriate for sim- 
ulations with small values of A. We have checked that 
results are consistent in such simulations by reducing S 
up to an order of magnitude. 

The deposited configurations are analyzed in search of 
arches. We first identify all mutually stable particles — 
which we define as directly connected — and then we find 
the arches as chains of connected particles. Two disks A 
and B arc mutually stable if A is the left supporting par- 
ticle of B and B is the right supporting particle of A, or 
viceversa. We measure the total number of arches, arch 
size distribution n(s), and the horizontal span distribu- 
tion of the arches n s (x). The latter is the probability 
density of finding an arch consisting of s disks with hor- 
izontal span between x and x + dx. The horizontal span 
(or lateral extension) is defined as the projection onto 
the horizontal axis of the segment that joins the centers 
of the right-end disk and the left-end disk in the arch. 



IV. RESULTS 

In Fig. 3 we show the area fraction (or packing frac- 
tion) cj) occupied by the disks as a function of the tapping 
amplitude. We also plot in Fig. 3 (z) SU pport as a function 
of A. We can identify four parts in Fig. 3(a). For small 
values of A (up to about 1.1), the area fraction falls very 
slowly from just above 0.84 to around 0.83. In the range 
1.1 < A < 1.15, simulations present large fluctuations 
but a sharp drop is observed; with 4> ranging between 
0.76 and 0.83. Then, in the range 1.15 < A < 1.5 there 
is a mild decrease in 0. Finally, for A > 1.5, we can 
observe in Fig. 3(a) a very slow rise in cj). 

The sharp change in </> around A = 1.13 has been de- 
scribed previously as an order-disorder transition due to 
the crystallization of the disks into a triangular lattice. 
This transition has been located around 4> = 0.8 by ex- 
periments |19| . which agrees with our simulations. The 
maximum value of <j> in the crystalline regime depends 
strongly on the ratio between the box width and the di- 



ameter of the disks. For integer ratios, the crystal can 
achieve the smallest lattice constant, hence the highest 
packing fraction. We do not commensurate our box to 
the particle size since this situation is difficult to achieve 
in experiments anyway. 

The slow increase in <p for very large tapping ampli- 
tudes is consistent with experiments carried out using 
non-circular particles "deposited" by a conveyor belt in 
a U-shaped container In this regime, we expect to 
slowly approach the value of sequentially deposited disks 
(i.e., (j) — 0-82 |21[). since a very strong tap should sepa- 
rate out the disks to such degree that they will fall back 
without many multi-particle collisions. 

In Fig. 4 we show two examples of packings with 250 
particles corresponding to the disordered [Fig. 4(a)] and 
the ordered [Fig. 4(b)] regimes. Arches are indicated 
with segments that join mutually stable particles. It is 
important to notice here that in our simulations particles 
that reach the base "stick" to their positions and are not 
pushed aside by further colliding disks during deposition. 
This prevents the formation of completely regular crys- 
tals since the first layer is already a randomly deposited 
layer. In a proper granular dynamic simulation one ob- 
tains more ordered structures j20j . Our ordered regime is 
better described as a layering of the system rather than 
a crystallization. This is indeed a limitation of the model 
since most 2D granular beds tend to present crystal-like 
order in real life. 

It is interesting to see whether in the transition region 
(1.1 < A < 1.15) the system presents a phase separation 
type behavior. We have plotted the profile of (j) m the 
vertical (y) direction for some tapping amplitudes in Fig. 
5. In the disordered regime, with exception of several 
bottom layers, the system shows a constant 4>(y) (with 
values below 0.8). In the ordered regime, the system has 
also a constant 4>(y) (with values around 0.84). Within 
the transition region, however, the profile shows a mono- 
tonic decrease from ordered-like area fractions down to 
disordered- like densities. We have found no sing of a step- 
wise profile which would indicate a coexistence between 
an ordered phase and a disordered phase. 

The most striking results from our simulations is the 
behavior of (.z)support as a function of A [see Fig. 3(b)]. 
Here, we can identify only three distinct parts. For small 
A (within the ordered region), (z) SU pport grows with A. 
This implies that larger tapping amplitudes "destroy" 
mutual contacts in the ordered phase. Interestingly, at 
the order-disorder transition, (z) support presents a sharp 
drop. This may be attributed to the appearance of many 
arches in the system that lead to the loss of order. Fi- 
nally, within the disordered region, the system increases 
its coordination linearly with A. This last observation 
is consistent with experiments in 2D of particles "de- 
posited" by a conveyor belt 0]. Clearly, arches are al- 
ways present in the system since (z) SU pport < 4 in all 
cases. The striking feature is the non-monotonic behav- 
ior of (z)support- Starting from small tapping amplitudes 
one can remove mutual contacts (and hence arches) by 
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increasing A. However, at the order-disorder transition 
a small increase in A leads to the creation of many new 
mutual contacts. Once in the new disordered regime, one 
can again remove mutual contacts by a further increase 
in A. Unfortunately, we are unable to support these find- 
ings with experimental evidence (apart from the behavior 
in the disordered regime that agrees with Blumenfeld et 
al. |l5j p since coordination numbers are rarely obtained 
in experiments. It is worth pointing out that simulations 
of 3D systems show a monotonic decrease (not increase!) 
in (z) as A is increased (see for example Ref. 

Form Eq. 0]we know that the number of arches in the 
system is inverse to (z) support , this is shown in Fig. 6 
where we plot the total number of arches per particle. 
Of course, we see again a sudden change at the order- 
disorder transition. Arches are more rare in the ordered 
phase just before the transition to the disordered phase. 
As we mentioned, in both regimes (ordered and disor- 
dered), arches are "destroyed" by increasing the tapping 
amplitude. 

It is commonly said that arches are responsible for 
voids in a granular pack which, in turn, diminish <fi. How- 
ever, from our results we can conclude that, in 2D, this 
is only true within the order-disorder transition region 
(1.1 < A < 1.15) and for very large tapping amplitudes 
(A > 1.5). In both cases a decrease in the number of 
arches correspond to an increase in the packing fraction, 
and viceversa. In contrast, for the rest of the packings 
(i.e., for A < 1.1 and 1.15 < A < 1.5), a decrease in 
the number of arches coincides with a decrease in pack- 
ing fraction. This result, may seem rather counterintu- 
itive; however, is not just the number of arches but their 
geometrical properties that determine the total volume 
of the voids left beneath them. Let us just remind our- 
selves that random packings of disks containing no arches 
at all (i.e., (z) S upport = 4) have (f) ~ 0.82 |2l|, which 
is lower than our arch-containing ordered packings and 
higher that our arch-containing disordered packings. 

In the inset of Fig. 6 we show the arch size distribution 
n(s) for a single tapping amplitude. This distribution is 
rather insensitive to A. However, very small changes in 
n(s) promote significant changes in (z) S upport- The arch 
size distribution can be fitted very well to a parabola in 
a semilogarithmic plot. However, this fitting can only be 
extended up to s = 7 or 8, since we find no arches with 
more than 8 disks in our simulations. Presumably, sim- 
ulations with a larger number of particles and a wider 
container will occasionally present larger arches. From 
Fig. 6, we can see that the simple exponential model we 
presented above is not suitable to describe arch distribu- 
tions in 2D packings. 

We have also analyzed the horizontal span of the 
arches. In Fig. 7 we represent n s (x) for s = 2 (a), 3 (b) 
and 4 (c), and the horizontal span n(x) averaged over all 
arches [Fig. 7(d)]. We include results from two simula- 
tions just before and after the order-disorder transition 
along with the theoretical model presented by To et al. 
based on a restricted random walk. In the case of 



n(x) we have no theoretical prediction since a theoretical 
model for n(s) is needed to weight the contributions of 
each size of arch. 

In the ordered regime, we can appreciate a clear dis- 
cretization of the arch extensions. This is easily under- 
stood since in the ordered state particles — even those 
forming arches — are organized in layers. Arches consist- 
ing of two disks, for example, can be formed by two disks 
at the same layer [corresponding to the peak near x = 1.0 
in Fig. 7(a)], or by one disk in one layer and another 
disk in the next layer displaced half " lattice constant" in 
x (corresponding to the peak near x — 0.5). This argu- 
ment can be extended to larger arches. See Fig. 4(b) to 
appreciate the form of the arches in the ordered regime. 

In the disordered regime, arches have no distinct lat- 
eral extensions but a continuous distribution. This is in 
agreement with the simple model by To et al. which was 
actually designed to represent arches at the outlet of a 
hopper. However, our arches tend to be more extended 
than those from the restricted random walk model. In 
particular, for two-disk- arches, we have no incidence of 
zero lateral span (corresponding to one disk on top of the 
other) in contrast with the model. 

Finally, in order to check the ability of the simulation 
model to reproduce some features typical of granular ma- 
terials, we have also performed an annealed tapping on 
the system. We have increased and decreased A progres- 
sively (from 1.0 to 2.0) five times and have plotted and 
(z)support as a function of A averaging all A-increasing 
cycles and all A-decreasing cycles separately. The very 
first A-increasing cycle is kept aside. Each increase from 
A = 1.0 to A = 2.0 takes 2000 taps. The results in Fig. 
8 correspond to averages over 5 independent simulation 
using 500 particles. For reference, we include the steady 
state curves obtained at constant tapping amplitude from 
Fig. 3. The inset shows results where the ramp rate is 
increased so that each cycle takes only 200 taps instead 
of 2000. 

According to Nowak et al. Q, after the very first in- 
creasing cycle, a granular bed should enter a reversible 
regime where further cycles of the tapping amplitude fol- 
low the same curve in the <p-A plot. However, Mehta 
and Barker (2^] found in their 3D simulations a hystere- 
sis loop and suggest that the area of the loop should 
increase as the ramp rate is increased. We also find hys- 
teresis in our 2D simulations. However, increasing the 
ramp rate does not clearly enlarge the hysteresis loop. 
The most clear change due to the increase in ramp rate 
is the difficulty that the system finds in reaching the or- 
dered regime. This is in agreement with the experiments 
by Nowak et al. where an increase in the ramp rate leads 
to lower densities in the small- A limit. 



V. CONCLUSIONS 

We have shown through a pseudo dynamic simulation 
that arches are ubiquitous in a 2D granular packing. Very 
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light tapping of the granular sample promotes layering of 
the particles. The transition from large-amplitude tap- 
ping (disordered packings) to small-amplitude tapping 
(ordered packings) is very sharp and occurs without a 
phase separation mechanism. Interestingly, we found 
that in a wide range of tapping amplitudes an increase 
in the number of arches coincides with an increase of 
the density of the packing. With exception of the order- 
disorder transition region, the mean coordination number 
is always increased by increasing the tapping amplitude; 



in contrast to what is found in 3D packings. 
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Figure captions 

FIG. 1. Mutual stability phase diagram for 2D pack- 
ings. Here, pi/N and P2/N are the fractions of particles 
that present 1 and 2 mutually stable supporting contacts, 



respectively. The state point of any packing must lie in 
the triangle defined by (0, 0), (0, 1) and (1, 0). The dotted 
lines correspond to states of equal (z) suppor t (the num- 
bers indicate the corresponding values). The solid curve 
represents all the states that can be obtained through the 
exponential model (see text). Open circles correspond to 
a representative sample of the packings generated in our 
simulations. 

FIG. 2. Packing fraction (a) and (z) suppc ,rt (b) as a 
function of the parameter S of the simulation algorithm 
for A — 1.1. Solid lines are only to guide the eye. 

FIG. 3. Packing fraction (a) and (z) SU p P ort (b) as a 
function of the tapping amplitude A. 

FIG. 4. Examples of the packings, (a) 250 particles 
tapped with A = 1.3 which corresponds to the disordered 
regime, (b) 250 particles tapped with A = 1.1 which 
corresponds to the ordered regime. Arches are indicated 
by segments joining mutually stable particles. 

FIG. 5. Vertical profile of the packing fraction for var- 
ious values of A: squares (1.1), circles (1.13), up trian- 
gles (1.135), stars (1-14), diamonds (1.15), down triangles 
(1.25), and pentagons (2.0). 

FIG. 6. Number of arches per unit particle as a func- 
tion of the tapping amplitude. The solid line is only a 
guide to the eye. The inset shows a semilog plot of the 
arch size distribution n(s). Symbols correspond to our 
simulations for A = 1.3 and the dotted line to a fit with 
a second order polynomial. 

FIG. 7. Distribution of the horizontal span n s (x) of 
the arches: (a) arches with two disks, (b) arches with 
three disks, and (c) arches with four disks. Part (d) corre- 
sponds to the horizontal span distribution over all arches. 
The dotted line corresponds to A = 1.2. The dashed line 
corresponds to A = 1.1. Solid line corresponds to the 
restricted random walk model (see Ref. Q). 

FIG. 8. Packing fraction (a) and (z) suppor t (b) as a 
function of the tapping amplitude A along an annealed 
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tapping. Tapping amplitude is increased and decreased 
several cycles. Each increasing (decreasing) cycle takes 
2000 taps. As a reference, we show the steady state ob- 
tained through constant tapping amplitude with symbols 
and solid thick line (see Fig. 3). The solid thin line corre- 
sponds to the initial increase of A. The dashed line corre- 



sponds to the decreasing phase and the dotted line to the 
increasing phase. All increasing cycles (apart from the 
initial increase) are averaged to produce a single curve. 
The same is done with the decreasing cycles. The insets 
show results where each increasing (decreasing) cycle is 
performed in only 200 taps. 



